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Abstract 

This chapter reviews numerical simulations of quantum field 
theories based on stochastic quantization and the Langevin equa- 
tion. The topics discussed include renormalization of finite step- 
size algorithms, Fourier acceleration, and the relation of the 
Langevin equation to hybrid stochastic algorithms and hybrid 
Monte Carlo. 



Invited chapter to appear in the special supplement "Stochastic Quan- 
tization" of Progress of Theoretical Physics. 



1 Introduction 



A great challenge in particle physics is the numerical solution of quantum 
field theories. Analogous problems appear in condensed matter physics as 
well. Stochastic quantization is extremely useful here, because it provides a 
direct path from the formal quantization of systems with many degrees of 
freedom to practical algorithms for numerical work. This chapter reviews 
this connection, leading up to the most popular algorithms for quantum 
chromodynamics. 

Quantum field theories involve an infinite number of degrees of freedom. 
This is the origin of ultraviolet divergencies. To perform any sensible calcu- 
lations, the number of degrees of freedom must be regulated. This can be 
done at the level of stochastic differential equations |]^, but for numerical 
computations one uses lattice field theory instead. The fields of continuous 
space-time are replaced by aggregate, or "block" fields, on the sites or links 
of a lattice For a recent review, see ref. If the lattice has a finite 
volume, lattice field theory is a quantum mechanical system with a large 
but finite number of degrees of freedom. 

One way to think of numerical simulations of lattice field theory is as 
Monte Carlo integration of the functional integral. An equivalent way is to 
imagine integrating the Langevin equations of stochastic quantization. The 
different viewpoints have led to different algorithms; with the latter approach 
leading to algorithms based on stochastic difference equations. As with de- 
terministic difference equations the general idea is to devise finite-difference 
approximations to the differential equations. However, the random noise 
affects the analysis of step-size errors in several ways. 

A central focus of this chapter is the dynamics of Langevin simulations. 
The term "dynamics" does not mean the dynamics of, say, QCD, but the 
behavior of the numerical algorithms in simulation, or CPU, time. In this 
language, the numerical algorithm is a dynamical system, whose static be- 
havior is the field theory (e.g. QCD) under study. To solve quantum field 
theories numerically, it is also important to analyze the dynamics of the al- 
gorithms, because fast dynamics obtain the solution in less computer time. 

This chapter is organized as follows. Sects. |2| and ^ introduce discrete 
Langevin and Fokker-Planck equations for scalar field theory. The analysis 
is the analog for stochastic differential equations of the analysis of step- 
size errors for deterministic differential equations. For local field theories 
renormalization plays an interesting role Q and the modifications needed 
for discrete Langevin equations are presented in sect. ^. It is argued that 
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the step-size errors of sects. ^ and |^ do not propagate to physical quanti- 
ties, when renormahzation is taken into account. For particle physics ap- 
plications non-Abelian gauge theories are the most important systems; they 
are treated briefly in sect. H. The analysis of these sections is extended 
to higher order integration schemes, such as Runge-Kutta in sect. |6[ The 
so-called hybrid stochastic algorithm is treated here. Sect. § explains algo- 
rithms for fermions. Sect. ^ discusses several schemes for accelerating the 
Langevin dynamics, i.e. for generating statistically independent lattice fields 
more quickly. Unfortunately, renormahzation does not eliminate step-size 
errors in algorithms for fermions or with acceleration, but sect. |9| presents a 
technique to make the algorithms exact. 

Langevin simulations for complex actions are covered in another chapter 

§■ 



2 Discrete Langevin Equations 

Consider scalar lattice field theory on a d-dimensional hypercubic lattice 
with spacing a. The standard action is 

S = a'^^(-i0,A(/>, + (2.1) 

X 

where Acpx := J^f^i^Px+fj. + (px-fi — 20a: )/«^; and a typical potential is 

= imV+5</''- (2.2) 

The subscripts x, etc, denote space-time coordinates. 

Stochastic quantization introduces a "time" parameter t. The fields 
evolve in f by a Langevin equation such as 

Mt) = -r]iit)-ViSit), (2.3) 

where the dot denotes a derivative with respect to t, the functional derivative 

and i is a multi-index denoting x and any internal indices. The noise is 
Gaussian: 

iVxit)) = 0, {r]x{t)7]y{u)) = 2 6{t- u)a-%y. (2.5) 
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The lattice spacing has been retained in eqs. (2.4) and (^]^) to determine 
the dimension of t in eq. ( p.3| ) . A scalar field has dimension [(j)] = ^{d — 2). 
Hence, [ViS] = i(d + 2) = [r/], whence [t] = -2. 
It can be shown that the t average 



0{(j)):= lim r dtO{(t){t)) (2.6) 



reproduces quantum mechanical vacuum expectation values, i.e. 



0{ct>) = ^ jmO{ct>)e-^. (2.7) 



In the lingo of numerical simulations, eq. ( |2.3| ) generates field configura- 
tions according to the distribution . 

In a numerical simulation t corresponds to computer time, but the com- 
puter can only evolve the fields in discrete steps e. A discrete approximation 
to the Langevin requires a finite-difference prescription for the derivative and 
a regularization of the Dirac (^-function. The Euler scheme is the simplest: 

</> ^ i (</>(* + e) - (/.(t)) (2.8) 

in eq. (|2.3D, and 

5{t -u)^ Utu (2.9) 

in eq. ( p.Sj ). Higher order schemes will be discussed in sect. ^. Writing 
X = t/e and r/^'^^ = r/(t/e)/y^ the Euler update becomes 

#^^^=#)-/f\ (2.10) 

where 

/f ) = V^r,f ) + eV.5W (2.11) 

and 

{vi':li'i)=2 6'^a-%y6at. (2.12) 

It is instructive to study the Euler algorithm in free field theory, where it 
can be solved exactly. Neglecting internal indices and Fourier-transforming 
to momentum space 

# = (l - ^^^(P))%f + V-ej: (l - eu'ip))'-'-' (2.13) 
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where uj'^{p) = + m? , p = 2sm^p. Now consider the correlations in 
Langevin time, which express the speed of convergence and de-correlation 
of the algorithm. From eq. ( |2.13| ) one finds 

4r'^4t^> = (i-^-^bi))''(4^^4t^> (2-14) 

where (•) denotes an average over the noise. Eq. ( 2.14| ) reveals many details 
of the performance of the Euler algorithm. First, taking the limit e — > 0, 
K — > oo with t = £K fixed, one sees that the convergence and de-correlation 
of cj)p occurs in time tc{p) = uj~'^{p). Second, for finite e one sees that the 
algorithm is stable only if |1 — ew^(p)| < 1 for all momenta. The most 
restrictive mode is the one with the largest momentum cj^j^^ — ^d/a? + m^. 
Once e is fixed by this ultraviolet mode, the infrared modes containing the 
interesting physics de-correlate in 

steps of the algorithm. An important characteristic of simulation algorithms 
is the critical dynamical exponent z, defined by N^. oc as a ^ 0, because 
the number of steps of the algorithm needed to completely de-correlate (p 
is determined by the largest N^- Typically z > 0, which means that more 
and more computation is needed to simulate field theories as the continuum 
limit is approached. This undesirable behavior is called critical slowing 
down. From eq. ( 2.15 ) one sees that z = 2 for the Euler algorithm. It will 
become clear in sects. ^, ^, and § that critical slowing down is closely related 
to the physical dimension of simulation time. Except for over-relaxation, cf. 
sect. |8.1| , the algorithms considered here all obey z = —[t] (for free field 
theory) . 



3 Discrete Fokker- Planck Equations 

One must now check that the probability distribution is correct as e — > 0. It 
is also useful to work out the 0(e) corrections to the distribution and develop 
a formalism for checking that higher-order schemes are indeed higher order. 
To simplify notation this section's equations are in lattice units, a = 1. 

Let the probability distribution at A be V^^'' [0]. The update of eq. ( |2.10| ) 
changes it to 

p(^+i)[</>]= j[d(t>'] /n5(0,-<A:+/.))^^'^[<^']- (3.1) 
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For small e the 5- functions can be expanded in powers of the drift force /j. 
Integrating over cj)' yields the Kramers-Moyal expansion 

oo 

= p(A)[^] + ^ _v,^ . . . V.„ • • • /.JT'Wm) . (3.2) 

n=l 

Eq. ( [3.2| ) gives a (functional) differential equation for the equilibrium distri- 
bution. Working to second order in e 



{fJ,) = 2e6ij+e^SiSj, 

-2 



(3.3) 



ififjfk) = 2e {SijSk + SjkSi + SkiSj) , 
ififjfkfi) = {^ij^ki + ^ik^ji + ^ii^jk) , 

using the abbreviation Si = ViS. 

To first order in e one obtains the Fokker-Planck equation 

V = Vi[{Si + Vi)V] (3.4) 

A change of variables V = e~^^^'^ brings eq. ( p.4| ) into the form 

i> = -{Si- V^) {S^ + Vi) ^ =: -n^. (3.5) 

In the space of (j) configurations, TC is self-adjoint and positive semi-definite. 
The unique zero-mode of Ti. is e~'^/^. Hence, generic initial conditions con- 
verge to the equilibrium solution ^ oc e~'^/^, or "P oc e^^ . 

Now let us analyze the next order in e. We are primarily interested in 
the equilibrium distribution, the solution to 

= Vi{S, + Vi)V + e{'^ViVjiSiSjV) + ViV\SiV) + iV^V^p} . (3.6) 

To simplify eq. (^) one repeatedly uses ViV = —SiP + 0(e) inside the 
braces, obtaining 

Q = Vi[{Si + Vi)Vl (3.7) 

where 

^ = ^ + 7 E i'^S,, - S]) . (3.8) 



4 
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Thus, the equilibrium distribution is P oc e , and S is called the equilib- 
rium action. 
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In field theory the detailed form of the action is not the whole story, 
because of renormalization. Indeed, the terms in the equilibrium action 
proportional to e are just those appearing in improved lattice actions. For 
example, changing field variables 

(Pi = ^^ + \eSi (3.9) 

changes the measure to 

[d(j)\^[d^]e^'^i^"'^ (3.10) 

and the action to 

s[^]^m+-jisi (3.11) 

i 

up to 0(e). Combining these two changes and writing = 
the probability distribution V oc e""^, where 

S = S + \Y.Sn- (3.12) 

i 

For polynomial actions this form modifies the bare couplings but induces no 
new terms. 

In Euclidean lattice field theory the interesting and accessible quantities 
are spectra and matrix elements. Since eq. ( |3.9| ) is a change a variables, 
it does not change the spectrum. Moreover, any multi-linear observable is 
correct up to second order: 

j [d$\ 0{4>) e-^ = J [#] O(0) e-^ + 0(^2), (3.13) 

since 

J [deb] Oj, = - 1 [#] V, (O,. e-^) = 0. (3.14) 

Hence, a trivial shift in bare couplings and a change of variables suffices to 
remove 0(e) terms from most interesting observables. 

4 The Detailed-Balance Universality Class 

The results of the previous section suggest that the non-zero step-size correc- 
tions do not affect physical predictions. Zinn- Justin has proven a theorem 
demonstrating the formal renormalizability of stochastic quantization |^ ^ . 
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The key ingredients of the proof are power counting, a BRST invariance, 
and a supersymmetry. This section examines how to adapt the arguments 
to analyze the non-zero step-size corrections to ah orders. One finds that 
BRST invariance still holds, which guarantees that the time discretization 
introduces only irrelevant operators. Then the critical phenomena of Euler 
and related algorithms is universal. But the supersymmetry does not hold, 
so the Fokker-Planck equation is not necessarily integrable. However, the 
pattern of supersymmetry violation is simple, and it is likely restored in the 
continuum limit. If so, the universality class includes the continuous-time 
Langevin equation. 

The mathematical steps are a direct transcription of ref. Q. Let us start 
with a general form of the dynamics, 

<,i = KM\- (4.1) 
The noise has a slightly different normalization than before 

«t> = 0, «t vl^) = e-Htu a-H.y (4.2) 

and the unit of time has been changed so that the step-size here is related 
to that of sects. ^ and ^ by e = 2e. The lattice spacing has been restored 
here, because we shall treat discrete space-time and discrete Langevin time 
on a similar footing. 

For convenience, let us introduce some notation. Greek letters will be 
used as a multi-index for space-time, internal indices, and Langevin time, 
viz. a = {x,a,t). For discrete time one must distinguish forward, backward 
and symmetric time difference operators, 

dt<t> = ±-^mt±e)-(t>{t)], (4.3) 

and 

di'U = Y,i^{t + e)-cP{t-e)]. (4.4) 

Note that dj-^"^ = ^{d^ + d^). Dot products, matrix products, or repeated 
Greek letters imply summation over all elements of the multi-index and 
multiplication of the sum by ea*^, e.g. 

J- (/> = Ja(/)a = eo"' V Ja(/'a. (4.5) 



By analogy with eq. (2^), the symbol = e a d / d4>'^{t). Finally, func- 



tional measures in this section are, for example, [dv] = Ha dva- 
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The generating functional for (dynamic) correlators of (j) is 

Z[J] = j [d(t>\[di'] 5{v - ^[0]) det>l e■^■'^-5^■^ (4.6) 
where the matrix M is given by 

= ypJ'a- (4.7) 

One represents the 5-function and determinant as functional integrals: 

5{v -T)= j [dC] e^-^"-^^ (4.8) 



with the contour of C,a along (— ioo,zoo), and 

detA^ = j [dc][dc]e^^^, (4.9) 



where the ghosts (c) and anti-ghosts (c) anti-commute. Inserting eqs. ( |4.8| ) 
and (O) into eq. (4^) and performing the (Gaussian) integral over u yields 

Z[J]= J[d^][dc][dc][dC]e-^^'^''''^''^^+-^-'^, (4.10) 

where the dynamical action 

S[^,c,c,C] = -^,C-C + C-^-cMc. (4.11) 
There is a BRST transformation, 

6Ca = 0, 
^Cq, = ^Ccf) 
5Ca = 0, 



(4.12) 



where ^ is an infinitesimal Grassman parameter. The terms of the dynamical 
action S transform as follows: 

5(-K-C)=o, 

HC-J')=^CMc, (4.13) 

6 (-cMc) = {(Mc + Ca V-yV/sJ^a CpC^) . 
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The last term vanishes because V^V/j^Fq, is symmetric under /3 <-> 7 whereas 
cpc^ is anti-symmetric. The two terms C,Mc cancel. Hence, the dynamical 
action is BRST-invariant. 

For the study of numerical algorithms a useful class of dynamics is given 

by 

^i,t = {dt<t^x,t + \Q%Sk^ (4.14) 



and Ai defined through T by eq. ([4.71). A Fokker-Planck analysis as in 
sect. 1^ shows that this Langevin equation converges to the correct proba- 
bility distribution, if is positive definite and independent of 0, and if 
lim,^o(Q^ - Q^) = 0. But we shall now apply the functional formalism to 
see what conclusions can be drawn for non-zero e. 

In the special case J^i^t = d^4^i^t + \^iS (i.e. the Euler scheme) there 
is an additional approximate symmetry. For discrete time the additional 
transformation is 

6(t)a = Cai, 

(4.15) 

5Ca = 0, 

for infinitesimal Grassman The terms of the dynamical action transform 
under eq. ( 4.15| ) as follows: 

<5(-iC-C) = 2c9j^^Ce, 

6iC-J^) = iCMc + 2di'h-J')t (4.16) 

<5 i-cMc) = i-cMC + 2cMd["U% 

Substituting the expressions for and M. into eq. ( [4. 161) and collecting all 
terms 

5S = 2c, [s^pdl'Up - d['^S^] i. (4.17) 

In a formal limit of continuous Langevin time the dynamical action is invari- 
ant, because dtSa = Sapdttpfs, from Leibniz' rule. For interacting theories 
with discrete Langevin time it is not, but the residue is a "lattice artifact." 
We shall return to this point after a discussion of renormalization. 

Ref. 1^] shows how power counting and the BRST invariance restrict 
the structure of the counter-terms. For the Euler update, where the step- 
size has (momentum) dimension [e] = —2 the dynamical fields have the 
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dimensions [C] = ^(d + 2), [</>] = i(d — 2) (as expected), and [c] + [c] = d. 
Renormalizability means that counter-terms in S of dimension greater than 
d + 2 are not needed, and the BRST symmetry relates ^-(p terms to c-c 
ones. The argumentation can be translated into the language of lattice field 
theory as follows. Any dynamics J- will have a relevant part of dimension 
^(d+2) and any matrix ea'^M. in the ghost action will have a relevant part of 
dimension 2. The BRST symmetry implies that the relevant parts of J- and 
M are related by eq. ( |4.7] ). In other words, there is a whole universality class 
of Langevin algorithms with the same critical dynamics. This universality 
class includes more sophisticated discretizations of Langevin time, such as 
those in sect. |^. In particular, the physics does not depend on e, up to the 
stability requirements discussed in sect. |2|. 

This conclusion is, perhaps, more easily digested by the following heuris- 
tic consideration. Restoring the lattice spacing, the Langevin step-size is 
e = eo^, by dimensional analysis. (The dimensionless number typed into 
the computer is e.) Consider a sequence of simulations with fixed e but the 
bare couplings of the static system (the model being simulated) tuned to ap- 
proach the continuum limit (of space-time) . Since simulation time is marked 
of in steps of la?, it seems to approach its continuum limit too. Hence, it 
is reasonable to guess that the non-zero step-size algorithms belong to the 
detailed-balance universality class, because the continuous Langevin equa- 
tion obeys detailed balance. This universality class includes the Metropolis 
algorithm and other exact, local algorithms. 

To prove the conjecture one must verify the probability distribution at 
equilibrium. From the non-perturbative proofs that stochastic quantization 
converges to the correct probability distribution one realizes that the 
supersymmetry plays an essential role. Therefore, let us return to the ap- 
proximate symmetry in eqs. (4.15)-(4.17). Even for Langevin time discrete. 



it combines with the BRST transformation to form something like a super- 
algebra. The generators D and D (of 5 and 5) satisfy = 0, = 
and 

{DD + DD)^ = 2d\'^^, (4.18) 

where <I> is 0, c, c, or C,. The right-hand-side is a discrete time-translation 
operator. Since 5S is a lattice artifact, it should be possible to adapt the 
approach of ref. Q to the supersymmetry of eqs. ( |4.12 ), ( 4.15 ) and ( 4.181) . 



(Ref. in proved that supersymmetry could be restored in the d = 2 Wess- 
Zumino model. The supersymmetry considered here is even simpler.) As- 
suming this strategy succeeds, the renormalized continuum limit of the dy- 
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namical theory is supersymmetric, and is the same as with detailed balance. 



5 Non-Abelian Spin and Gauge Systems 

In particle physics the most interesting field theories are non-Abelian gauge 
theories and chiral models. On the lattice the fundamental variables are 
Lie group elements defined on links (gauge theories) or sites (chiral models) 
of the lattice. The analysis of the previous sections can be adapted to 
non-Abelian theories. The crucial ingredient is to define differentiation in 
the group manifold in a way consistent with partial integration over Haar 
measure. With such a definition of Vj, eq. (3.2) still holds. 



We shall concentrate on unitary groups and use the following conven- 
tions: The anti-Hermitian generators T"" are normalized by Tr(r"r^) = 
— i(5"*. The structure constants are given by the commutation relations 
[r'^,r^] = -f^bcrpc^ g^^j^ parameters and write ixi = uo'^T'^. The 

derivative is defined by [0] 

fie^U) = f{U) + u'^WiU) + 0{u;^), (5.1) 

where f{U) is any function of the unitary matrix U. The most useful ex- 
ample is 

VTviUV) = TriT^UV), VTriVU^) = - TriVU^T"), (5.2) 

where V is independent of U. The derivatives do not commute (the Lie- 
group manifold is curved), [V",V^] = —/"-^^V^. The commutation relation 
is especially easy to verify from eq. (|5.2D . 

In field theory one must keep track of a collection of unitary-group de- 
grees of freedom. The commutation relation reads 

[K,>.My,u]=-r'''^'S^y^u. (5.3) 

for gauge fields, and a similar expression without the labels fi, v for spin 
fields. It is convenient to introduce a multi-index i = (a, x, \x) for gauge 
fields and i = (a, x) for spin fields. In the following we shall concentrate on 
gauge fields. 

A Langevin update for non-Abelian fields is given by 

C/if;i) = e-/^%^^C/W. (5.4) 
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The Euler drift force is 

fi = ^/er|i + eViS, 



(5.5) 



where r]i = f]'^ ^ are Gaussian random numbers with dispersion 2. 

The equihbrium action can be worked out just as in sect. |3|, taking care 
that the Vj now longer commute. This leads to an new 0(e) "correction" 
to the equilibrium action, which now reads 

m = (l + ^) S[U] + \Y. {^"^Pp] - {V,S[U]f] , (5.6) 

j 

where Ca is the Casimir invariant of the adjoint representation {Ca = N for 
SU(A^)). The S'j term can again be absorbed into a change of variables, and 
for a simple plaquette action the remaining 0(e) terms are absorbed into 
the bare couplings |T^. For example, the shift in (5 of the Wilson action is 
(5 ^ j3[\ + e(CA/12 — Cp)], taking the shift from eq. ( |5.6| ) and the change of 
variables into account Wilson loops are multi-linear, so that, after the 
change of variables, their expectations are correct up to 0(e). Alternatively, 
one can correct a Wilson loop by a factor of 1 + eCp/^ per link. Note that 
this correction factor drops out of a Creutz ratio. 

The analysis of sect. Q can also be extended to (pure) gauge theories. 
Once again, because [e] = —2, one expects that dynamical systems with 
different values of e belong to the same universality class, i.e. that renor- 
malization washes out non-zero step-size effects. 

6 Higher Order and "Hybrid" Algorithms 

For the systems considered so far, the effects of discrete Langevin time are 
analogous to lattice artifacts. However, it is still sometimes desirable to 
investigate discretizations that suppress them. For example, to find a non- 
trivial fixed point one must investigate the phase structure of the lattice 
theory; this may require more precise control over parameter-space than 
what the Euler update would allow. Also, in the following sections we 
shall investigate dynamics for which the renormalization theorem does not 
apply. For example, to eliminate critical slowing down, it is necessary to 
introduce dynamics with dimensionless time and, hence, different power 
counting. Furthermore, fermionic systems are almost always treated by a 
"pseudo-fermion" fields whose interactions with scalar or gauge fields is non- 
local and non-polynomial. 
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This section considers algorithms that are still approximate, but have 
smaller 0(e) effects. Exact algorithms are considered in sect. |9[ 

First, we shall consider the Runge-Kutta algorithm [|l2|, |l^. The new 
configuration is obtained from the old one by 

U^T^ = e-f-'-U^i, (6.1) 

with 

/. = V^V^ + + \eCA) {Sf''^ + S[^^^'^^) (6.2) 

where S^^'^^^'^^ denotes the action evaluated using the tentative update 

f;(A+l/2)^^-/.„^(A) (g_3) 

where / is an Euler update with the same noise as in /. Expanding S^'^^^^'^^ 
in powers of ^/e and working out the changes to eqs. (|3.3| ) and (|3.6|) , one 
finds the equilibrium action coincides with the desired action up to terms of 
O(e^). 

Another update with O(e^) accuracy is obtained by eq. (6T) with 

h = V^Vi + ^(1 + i-,eCA)S, - ie3/2 ^ Sij7]j (6.4) 



and no tentative update. For scalar or pure gauge theories eq. ( p7i| ) has 
no advantage over the Runge-Kutta scheme, which is easier to implement. 
However, the generalization when fermions are coupled in saves an expensive 
matrix inversion at each step. 

The standard hybrid stochastic algorithm can also be considered as an 
improved discretization of the Langevin equation. Consider the following 
update steps for scalar field theory: 



^(A,l/2)^^(A,0)_i^^(A,0) 
.(A,1)^.(A,0)^ (A,l/2) 



(6.5) 



where Gaussian noise with unit dispersion, followed by a molecular 

dynamics p3|] trajectory of — 1 steps of 



(A,n+l/2) _ (A,n~l/2) . ^(A,n) 
,(A,n+l)^ ,(A,n) ^^^(A,n+l/2) 



(6.6) 
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Finally, the new configuration is given by 



, (A+1,0) _ ,{A,iV) 



(6.7) 



Eq. ( |6.6D is the leap-frog scheme for integrating Hamilton's equations for 
Hamiltonian H = \ Yl,i '^i + ^[4'] through the (vr, (p) phase space; eliminat- 
ing vr it is equivalent to the Verlet scheme for integrating a second-order 
ordinary differential equation. Eq. ( |6.5| ) is the occasional "refreshment" of 
the velocities needed to average over momentum. 

For = 1 the hybrid update collapses to the Euler discretization of the 
Langevin equation. The is the basis of the statement that "Langevin is a 
special case of hybrid." For N = 2 the hybrid update is similar in structure 
to the Runge-Kutta update, but instead of O(e^) accuracy the equilibrium 
action is given by eq. ( |3.8| ) with e = For arbitrary the structure is 
similar to higher-order Runge-Kutta schemes, but the equilibrium action is 
the same for all A^. 

Again it is instructive to analyze the performance of these algorithms in 
free lattice field theory. In momentum space the leap-frog iteration can be 
worked out 



where sin^ and cos at are polynomial approximations of sine and cosine: 

siuM = E (-1)"7^^ n 1 - (6-9) 



n=0 ^ ' m=0 

and 



N a2n n-1 2 

cosM^) = 5:(-ir— - ni-]^- (6-10) 



n=0 ^ m=0 

The auto-correlation function is then 



= co^UN5u{pi)) , (6.11) 

where (•) denotes an average over the every noise vr^''*''^^ 

Let us consider two idealized limits. One is the "Langevin limit" 

fixed, K ^ oo, 5^0, t = \kN'^5'^ fixed, (6.12) 
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and the other is the "molecular dynamics limit" 

K fixed, iV ^ cx), 5^0, t = N5 fixed. (6.13) 

In the Langevin limit, cos'^{N5uj) e"*"^^ and the dynamics de-correlates 
as any Langevin dynamics. In particular, z = 2. In the molecular dynamics 
limit cos^(A^(5w) cos'^ (lot). If ujt small this limit differs little from the 
Langevin limit, and in particular the computation needed to de-correlate 
the slow modes is about as much as with the Euler update. If one tries 
to make the trajectory length longer, a problem arises because there is a 
spread of frequencies in field theory. For r > 27r/a;niax almost every choice 
of T coincides with a multiple of 2'Kn/uj{p) for some p and n [14|. This is a 



piece of bad luck, because such a mode never de-correlates. As the physical 
volume of the system increases the density of modes increases until r has 
nothing but bad luck at all. Under these circumstances it is difficult to 
define z sensibly. 

It is safe to say that an important attraction of the hybrid algorithm 
was the claim that it had z = 1 because the time parameter of the tra- 
jectories had (momentum) dimension [r] = [S\ = —1. However, numerical 
studies have shown that the optimal trajectory length is Topt ~ 7r/2a;niax5 in 
accord with the above remarks. Therefore, the fast, ultraviolet modes de- 
correlate in one trajectory, but the slow, infrared modes de-correlate as in 
usual Langevin dynamics. For the standard hybrid algorithm the short, fixed 
trajectory length chosen makes it nothing but an elaborate discretization of 
the Langevin equation, with step-size e' = ^r^. This step-size is larger 
than in the Euler algorithm (for equal step-size error), but nevertheless the 
stochastic process has dynamical critical exponent z = 2. When step-size 
errors matter, it is not clear which discretization of the Langevin equation 
is preferable, hybrid or Runge-Kutta, when all aspects of the computation 
are considered. 

The trajectory length can be increased to roughly 2ti/ LOmm if its length 
varied from trajectory to trajectory ||l^. Remarkably, this solution is an 
element of the original hybrid algorithm |jl^, in which the trajectory length 
was to be chosen with probability (1 — P5)^ . Then, although any given 
mode has bad luck occasionally, most of the trajectories de-correlate it. With 
variable trajectory length and the option to select any configuration 
for the ensemble, the stochastic process is no longer in the detailed-balance 
universality class. In particular, the proof of convergence must be modified 
]l6| , p!7[ and the formalism of sect. ^ does not apply. Nevertheless, the 
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equilibrium action is still given by eq. (|3.8|) (or eq. (|5.6D for gauge theories) 



with e = 1 17]. An individual harmonic oscillator has auto-correlation 



time Tc = 2/P, provided P <2uj |T^. For free field theory it is then easy to 
see that choosing P = 2u;mm de-correlates all modes in (molecular dynamics) 
time Tc = As in standard Langevin the maximum step-size is set by 

stability of the ultraviolet modes. Given this step-size, the number of sweeps 
needed to de-correlate the infrared modes is 

A^e = T 0^ — • (6-14) 
am 

The dynamical critical exponent has been reduced to z = 1, because for 
random trajectory lengths, the parameters can be chosen so that CPU time 
is molecular dynamics time. 



7 Including Fermions 

The previous sections considered systems with bosonic degrees of freedom 
only. This section treats algorithms for fermionic degrees of freedom. 

First note that the fermions' part of the action can always be formulated 
in a quadratic form. Then 

:#][#]e^*^[^l^ = detM[;7], (7.1) 

where U denotes gauge fields interacting with the fermions. For QCD, Af [[/] 
is a lattice version of ^ + m, so it is not real. The standard remedy is to 
introduce a second flavor of fermion and introduce a complex field |18|: 

detM[J7] detM[;7] = detMt[;7] detM[;7] = j [d^^Je-V, (7.2) 

where Spf = {M^ [U]M[U])-^ if , and 93 is often called the pseudo-fermion 
field. Since = 75M75 the pseudo-fermion action can be written 



5pf 



(^tM^-V, (7.3) 



where M5 = 75M. Below these details are less important the the form of 
the action in eq. (^^), so the subscript 5 will be dropped. 
The Euler update for the pseudo-fermion is [^, 11] 



f = - e^Mr^) V^f + V^C^ (7.4) 
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where i is a multi-index for space-time, color, and spin indices of ^p, and the 
Gaussian noise has dispersion (^^^) = 2. The drift force of a gauge field 
coupled to the fermion is augmented by a new term 

fi^fi- evptM-2(ViM~2)M-V- (7.5) 

Notice that the step-sizes E^p and e appearing in the fermion and gauge 
updates need not be the same. 

These dynamics suffer from a peculiar critical slowing down. The eigen- 
values of for free Wilson fermions are //(p) = sin^(pa)/a^ -|- {m + ^ap^)^. 
The fastest modes in eq. ( |7.4| ) are the low momentum modes, /U~^(0) = m~^, 
and their stability restricts the magnitude o^£^p. The auto-correlation length 
of the high momentum modes is long, = {m + 2d/a)'^, and consequently 
they de-correlate in Nc oc {am)~'^ sweeps. 

However, it is easy to eliminate this critical slowing entirely. Consider 
eq. (PI ) with Q = Q = M, i.e. 



V^f = (1 - e^)(^f ) + M,,^Ci, (7.6) 

and (^t^) = 2- e^. One can show from the Fokker-Planck equation (or 
BRST techniques!) that the equilibrium distribution of ip is correct to all 
orders in e^, with this modification in the noise. One can even set = 1, 



in which case the fermion field de-correlates immediately [11|. Then the 



gauge-field drift force is augmented by the bilinear noise term 

fi^fi- ee^Ae, (7.7) 

where A = M-\ViM-^)M-\ and (^t^ = i. 

The bilinear noise algorithm of eq. (^]^) can be extended to a higher- 
order scheme using variations of the Runge-Kutta technique ||2^, ^] . How- 
ever, there is a complication. Terms of the form (^^AC^^A'C) ™^ke it im- 
possible to integrate the Fokker-Planck equation. Unfortunately, the higher- 
order step-size errors in procedures that remove the non-integrable terms are 



proportional to the volume [22|. Hence, even though the error is formally 
higher-order, the step-size must be chosen to be smaller. On the other hand, 
numerical work ||2^ indicates that the non-integrable terms make only a 
small contribution to observables, and it is less harmful to leave them in. 

In combining the fermion updates into hybrid algorithms, several ap- 
proaches are possible. One can leave ip = fixed during the gauge- 
trajectory, update (f for fixed ^, or generate a new ^ at each step of the 
trajectory. For an analysis of these possibilities, see ref. [p^]. 
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8 Accelerating the Dynamics 



In free field tlieory the original hybrid algorithm ameliorates, but does not 
eliminate, critical slowing down. This section uses the Langevin equation to 
explore two other ways to attack the problem, over-relaxation and Fourier 
acceleration. In the former case the Langevin equation is used as a peda- 
gogical tool; most implementations rely on other algorithms. In the latter, 
however, stochastic difference equations are essential. 

The name of the game is to accelerate the dynamics of the slow modes 
and thereby reduce the critical dynamical exponent. It can be determined 
analytically in free field theory, but reliable determinations for strongly in- 
teracting systems are extremely difficult. For four-dimensional interacting 
systems, such as QCD, it has not yet proven feasible to quote z with sensible 
error estimates. 



8.1 Over-relaxation 

Let us first consider over-relaxation [^] . The original formulation and most 
practical implementations do not look much like the Langevin equation. It 
is, however, possible to re-cast it into this form |25]. Imagine two harmonic 
oscillators, i.e. action S = + \'^\<^\- A Langevin equation with the 

properties of over-relaxation is 



the Langevin dynamics couples the two modes together. More generally, 

4)% = —(cos 6 5ij + sin 6 f-ij)Sj + V cos 9 rji , (8-2) 

so the mode-coupling term drops out of the Fokker-Planck equation. (Here 
eij is the anti-symmetric tensor.) Hence, eq. (|8.1| ) generates configurations 
with the correct probability distribution, if cos 9 > 0. 





The eigenvalues of the matrix in eq. (8T) dictate convergence. They are 



= i(u;?+L^|) cos ± iy'(a;2 _ ^2)2 _ (^2+^2)2 sin2^ (8.3) 
Clearly, if 



2^2 
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the eigenvalues form a complex conjugate pair with Rei/± ^ iOii02- The off- 
diagonal coupling in eq. accelerates the slow mode at the expense of 
decelerating the fast mode. In free lattice field theory uj'^{p) = + m?, and 
a typical strategy couples a mode with momentum p to one with momentum 
P = (tt + pi, TT + •••)• Then = ^djc? — ■ The angle Q is chosen so 
that eq. ( ^.41) holds for all p. For example, choosing sin0 = 2d/{2d + a^m?) 
one finds T~^{p) = Rez^± = {m/a)\/4d + a^vn? independent of p. In typical 
implementations a dimensionless step-size is fixed. Hence, the de-correlation 
time measured in sweeps is A^c oc (om)~^. Like the original hybrid scheme, 
over-relaxation reduces the dynamical critical exponent to z = 1, but does 
not eliminate critical slowing down completely. 



8.2 Fourier Acceleration 

In free lattice field theory Langevin updating can be studied exactly in 
momentum space, cf. sect. |2[ This analysis shows why the usual dynamics 
have z = 2. It also suggests a remedy. In eq. ( ^.13 ), instead of taking the 



step-size independent of p, one could just as well take e{p) = e/uj {p). This 



is Fourier acceleration |2C, |lT|. The natural time step is now e, which is 
dimensionless, so (for free field theory) it is easy to see that all modes de- 
correlate on the same time scale, and that that time scale is independent of 
a. This would eliminate critical slowing down at the theoretical level. The 
more relevant standard of success is computation. Fortunately, the cost of 
the fast Fourier transform (FFT) algorithm increases only as V log V . 

For interacting theories the central question is the tolerable value of e. 
In position space e becomes non-local, e = eQ^y^ where 

2 f d'^p e'^P^^~y'> 
'^"^ " J (2^ p2 + m2- ^^-^^ 

To leading order in e the equilibrium action becomes |11] 

S = S + ^Y.Qh('^S,j-S,Sj). (8.6) 

These interactions can be made local by introducing a new field ( with ki- 
netic term ^({A + m?)( and cp-C interactions QSi and QSijCj with couplings 
proportional to e. With a local field theory, familiar techniques can be used 
to predict the form of C interactions on physical quantities. (It is difficult 
to determine their size, except that they are proportional to e.) Once the 
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form is known, the step-size errors can be eliminated by extrapolating — in 
essence one takes the continuum limit of Langevin time explicitly. 

The field theoretic analysis can proceed from eq. ( |8.6| ) , or one can use the 
formalism of sect. ^ Starting from the dynamical action in eq. ( 4.11 ), it is 



convenient to change variables ( ^ CQ~^ c i— > cQ~^, to make the theory 
local. The algebraic form of the BRST transformation does not change, so it 
can be used to derive identities relating different quantities, e.g. correlators 
with and without (^'s and ghosts. In particular, the C field has space-time 
interactions; it is the same field as in the previous paragraph. 

If predictions of the e dependence are not reliable enough to extrap- 
olate, generalizations of the Runge-Kutta method are available |l^, even 
when fermions are included |2^, However, while Runge-Kutta processes 
render the step-size errors O(e^), the remaining errors are too cumbersome 
to analyze conveniently. 

A complication arises for non-Abelian gauge theories. The eigenvalues of 
the covariant (lattice) Laplacian are approximately labeled by momentum 
only in a smooth gauge. Consequently, it is necessary to fix the gauge before 
applying Fourier acceleration, which makes implicitly field-dependent. 
An alternative, using covariant derivatives in the definition of Q^, makes it 
explicitly field-dependent. Either method changes the Fokker-Planck equa- 
tion at leading order. The first two moments in eq. (|3.3|) become 



(/,)=e-Q2.5,, 



•7) 



{fdj) = 2eQl. 
The Fokker-Planck equation is then 

V = ViQl [{Sj + Qj.'ViQli) + VjV 

which equilibrates to the wrong distribution. This must be repaired by 
replacing Sj with Sj — logQ^; in the drift force. In practice, the repair is 



implemented using a stochastic estimator, cf. ref. |27|. 



9 Exact Algorithms 

For systems with fermions, e.g. QCD, the non-zero step-size errors can alter 
physical results. Although the effects can be analyzed, the analysis is not 
especially straightforward, because the pseudo-fermions interact non-locally. 
For Fourier accelerated algorithms, the step-size errors are easier to analyze. 
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because the algorithm can be re-cast as a local theory. In both cases, how- 
ever, an exact algorithm is desirable. Even for models where renormalization 
is thought to wash out step-size effects, most people would prefer an exact 
algorithm, at least for psychological reasons. 

There is an exact algorithm, well-suited to QCD, called the hybrid Monte 



Carlo [^. It is based on the hybrid scheme discussed in sect, g, but configu- 



rations are accepted or rejected according to a Metropolis test. The secret is 
to apply the test to H[tt, (p] = ^vr^ -|- £'[(/>], rather than to S[(j)] alone. Starting 
with the configuration <J)^^'^\ the steps are as follows: 

1. Generate vr^'^''^^ from a Gaussian distribution e~^^ . 



2. Carry out eq. ( |6.5| ) and — 1 steps of eq. ( |6.6[ ). 

3. Bring ir up to the same (molecular dynamics) time as (p: 

^iX,N) ^ ^iX,N-l/2) _ ,^giX,N) ^ (9.1) 

4. Make the substitution = with probability min (l, e~^^^ , 
where := H^^'^^ - ij(^'O); otherwise = unchanged. 

When the process is iterated the configurations labeled have the de- 

sired probability distribution V = . As in sect. an improvement is to 
randomize the value of in step |2| . 

The hybrid Monte Carlo algorithm has become the algorithm of choice 
in numerical simulations of full QCD. As such it warrants a review of its 
own, but that is beyond the scope of a set of articles on stochastic quan- 
tization. For a review see ref. [^ ]. To give fair comparison to the other 
algorithms, however, let us work out the critical dynamical exponent (for 
free field theory). The number of trajectories needed to de-correlate the slow 
modes is A',- oc (roptWinin)"^ and the amount of computation in a trajectory 
is proportional to = Topt/(J. Therefore, the total amount of computation 
needed to de-correlate the slow modes is 

= N,N oc oc (9.2) 

since (Jwmax ^ 1 for stability of the molecular dynamics trajectory. For a 
fixed trajectories of length Topt ~ "^maxi one sees that A^c oc (am)~^, i.e. 
z = 2. Similarly, for trajectories of variable length with mean Topt ~ "^min' 
critical slowing down is less severe, and z = 1. 
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Hybrid Monte Carlo has an additional source of slowing down in the 
infinite volume limit with lattice spacing fixed. A leap-frog trajectory drifts 
off the energy sheh by an amount of size |0| AH = C^b'^ + iCV(5*^ + • • •. 
Consequently, one must reduce the step-size as (5 oc F"^/^, otherwise the 
Metropolis test in step [^rejects almost every trajectory. This infinite volume 
slowing down can be alleviated by higher-order integration schemes. Since 
the number of degrees of freedom increases in the continuum limit too, this 
characteristic could affect the values given for z in free field theory. 

For QCD hybrid Monte Carlo, with random trajectories and some of the 
other ideas proposed in ref. [31|, seems to be the algorithm of the near future. 
In practice the approach to the lattice-spacing and volume limits is restricted 
by computer memory. As a four-dimensional theory, QCD has enormous 
memory demands, so with only moderate critical slowing down and despite 
some infinite-volume slowing down, an appropriately tuned hybrid Monte 
Carlo should be adequate. 
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